Chapter 4 Community composition
4.1 Taxonomy overview
4.1.1 Stacked barplot
# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
filter(count > 0) #filter 0 counts
# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)
ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ type, scales="free") + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
labs(fill="Phylum",y = "Relative abundance",x="Samples")+
scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label
4.1.2 Phylum relative abundances
phylum_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
group_by(sample,phylum) %>%
summarise(relabun=sum(count))
phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| phylum | mean | sd |
|---|---|---|
| p__Bacteroidota | 0.380806816 | 0.202118047 |
| p__Bacillota_A | 0.253813952 | 0.163904042 |
| p__Bacillota | 0.118079729 | 0.146679823 |
| p__Pseudomonadota | 0.096901291 | 0.160581847 |
| p__Campylobacterota | 0.053050348 | 0.092998130 |
| p__Verrucomicrobiota | 0.029867960 | 0.073128363 |
| p__Desulfobacterota | 0.023580475 | 0.036482484 |
| p__Chlamydiota | 0.010572704 | 0.059690269 |
| p__Fusobacteriota | 0.010482132 | 0.028311738 |
| p__Cyanobacteriota | 0.009206465 | 0.016492133 |
| p__Bacillota_C | 0.004713164 | 0.006645703 |
| p__Spirochaetota | 0.004009680 | 0.012308028 |
| p__Bacillota_B | 0.002465465 | 0.004927667 |
| p__Actinomycetota | 0.001235741 | 0.006346852 |
| p__Elusimicrobiota | 0.001214079 | 0.006501752 |
phylum_arrange <- phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun)) %>%
arrange(-mean) %>%
select(phylum) %>%
pull()
phylum_summary %>%
filter(phylum %in% phylum_arrange) %>%
mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
geom_jitter(alpha=0.5) +
theme_minimal() +
theme(legend.position="none") +
labs(y="Phylum",x="Relative abundance")
4.2 Taxonomy boxplot
4.2.1 Family
family_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,family) %>%
summarise(relabun=sum(count))
family_summary %>%
group_by(family) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| family | mean | sd |
|---|---|---|
| f__Bacteroidaceae | 2.215146e-01 | 0.1392048812 |
| f__Lachnospiraceae | 1.406192e-01 | 0.1085410432 |
| f__Tannerellaceae | 1.028077e-01 | 0.0799387840 |
| f__Helicobacteraceae | 5.258990e-02 | 0.0925721039 |
| f__Mycoplasmoidaceae | 3.694491e-02 | 0.0756423660 |
| f__Erysipelotrichaceae | 3.502391e-02 | 0.0451216980 |
| f__UBA3700 | 3.418595e-02 | 0.0557485850 |
| f__Marinifilaceae | 2.774048e-02 | 0.0271999075 |
| f__ | 2.772765e-02 | 0.0838376915 |
| f__DTU072 | 2.700317e-02 | 0.0529029477 |
| f__Rikenellaceae | 2.696900e-02 | 0.0464477971 |
| f__Enterobacteriaceae | 2.691677e-02 | 0.0918386902 |
| f__Coprobacillaceae | 2.551773e-02 | 0.0892179821 |
| f__Desulfovibrionaceae | 2.358047e-02 | 0.0364824837 |
| f__Ruminococcaceae | 1.859938e-02 | 0.0429483670 |
| f__LL51 | 1.759437e-02 | 0.0687507595 |
| f__Rhizobiaceae | 1.596807e-02 | 0.0770633474 |
| f__UBA3830 | 1.578817e-02 | 0.0454016753 |
| f__Akkermansiaceae | 1.227359e-02 | 0.0316838689 |
| f__Chlamydiaceae | 1.057270e-02 | 0.0596902690 |
| f__Fusobacteriaceae | 1.048213e-02 | 0.0283117384 |
| f__CAG-239 | 9.141331e-03 | 0.0150992357 |
| f__Enterococcaceae | 8.420523e-03 | 0.0466580889 |
| f__Gastranaerophilaceae | 7.728770e-03 | 0.0144404459 |
| f__Oscillospiraceae | 6.643218e-03 | 0.0078105983 |
| f__UBA1997 | 6.378408e-03 | 0.0309668631 |
| f__Streptococcaceae | 6.366441e-03 | 0.0342174908 |
| f__UBA1242 | 4.673359e-03 | 0.0153730061 |
| f__Brevinemataceae | 4.009680e-03 | 0.0123080282 |
| f__Acutalibacteraceae | 3.374550e-03 | 0.0109560112 |
| f__UBA660 | 3.196511e-03 | 0.0139558713 |
| f__Clostridiaceae | 2.842205e-03 | 0.0171339212 |
| f__RUG11792 | 2.817729e-03 | 0.0251127757 |
| f__Peptococcaceae | 2.465465e-03 | 0.0049276669 |
| f__MGBC116941 | 2.160169e-03 | 0.0093827415 |
| f__Acidaminococcaceae | 1.943804e-03 | 0.0050321931 |
| f__CAG-508 | 1.807978e-03 | 0.0064175379 |
| f__Anaerovoracaceae | 1.569531e-03 | 0.0036471440 |
| f__Moraxellaceae | 1.485415e-03 | 0.0097446552 |
| f__RUG14156 | 1.477695e-03 | 0.0045847831 |
| f__Staphylococcaceae | 1.361709e-03 | 0.0050863070 |
| f__Elusimicrobiaceae | 1.214079e-03 | 0.0065017516 |
| f__CAG-288 | 9.490865e-04 | 0.0060198775 |
| f__Anaerotignaceae | 8.989745e-04 | 0.0040469504 |
| f__CALVMC01 | 7.516846e-04 | 0.0043609744 |
| f__Eggerthellaceae | 6.407882e-04 | 0.0021266467 |
| f__Massilibacillaceae | 6.235073e-04 | 0.0016350480 |
| f__Mycobacteriaceae | 5.949530e-04 | 0.0060400054 |
| f__UBA1820 | 4.736030e-04 | 0.0013057353 |
| f__Arcobacteraceae | 4.604457e-04 | 0.0050324221 |
| f__CAG-274 | 4.519746e-04 | 0.0022028477 |
| f__Burkholderiaceae_C | 3.699430e-04 | 0.0048092594 |
| f__Muribaculaceae | 3.617894e-04 | 0.0009776409 |
| f__UBA932 | 3.419549e-04 | 0.0011583178 |
| f__Hepatoplasmataceae | 2.989107e-04 | 0.0038858387 |
| f__Rhodobacteraceae | 2.959092e-04 | 0.0038468195 |
| f__Weeksellaceae | 2.771627e-04 | 0.0031476408 |
| f__Eubacteriaceae | 1.646822e-04 | 0.0006729067 |
| f__Sphingobacteriaceae | 1.505775e-04 | 0.0012460017 |
| f__Devosiaceae | 1.489995e-04 | 0.0015094318 |
| f__Pumilibacteraceae | 1.277418e-04 | 0.0007646754 |
| f__WRAU01 | 9.603359e-05 | 0.0012484367 |
| f__Peptostreptococcaceae | 2.287338e-05 | 0.0002973539 |
family_arrange <- family_summary %>%
group_by(family) %>%
summarise(mean=sum(relabun)) %>%
arrange(-mean) %>%
select(family) %>%
pull()
family_summary %>%
left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
filter(family %in% family_arrange[1:20]) %>%
mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
scale_color_manual(values=phylum_colors[-8]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")
4.2.2 Genus
genus_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,genus) %>%
summarise(relabun=sum(count)) %>%
filter(genus != "g__")
genus_summary %>%
group_by(genus) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| genus | mean | sd |
|---|---|---|
| g__Bacteroides | 1.352836e-01 | 0.0923254347 |
| g__Parabacteroides | 9.681305e-02 | 0.0801790006 |
| g__Phocaeicola | 6.935126e-02 | 0.0794715561 |
| g__Mycoplasmoides | 3.051289e-02 | 0.0753256189 |
| g__Helicobacter_J | 3.009004e-02 | 0.0594808174 |
| g__Odoribacter | 2.552307e-02 | 0.0267570773 |
| g__Roseburia | 2.384718e-02 | 0.0559439989 |
| g__NHYM01 | 2.249986e-02 | 0.0796951054 |
| g__Alistipes | 2.208351e-02 | 0.0283591693 |
| g__Coprobacillus | 2.013859e-02 | 0.0878851683 |
| g__Agrobacterium | 1.596807e-02 | 0.0770633474 |
| g__Akkermansia | 1.227359e-02 | 0.0316838689 |
| g__Fusobacterium_A | 1.038903e-02 | 0.0283170720 |
| g__Kineothrix | 8.801109e-03 | 0.0409037874 |
| g__Proteus | 8.657984e-03 | 0.0681831315 |
| g__Dielma | 8.578108e-03 | 0.0090836958 |
| g__CAG-95 | 8.109350e-03 | 0.0205019603 |
| g__JAAYNV01 | 7.994602e-03 | 0.0195404994 |
| g__UBA866 | 7.260325e-03 | 0.0293499698 |
| g__Desulfovibrio | 7.018155e-03 | 0.0211479806 |
| g__Enterococcus | 7.001300e-03 | 0.0456600395 |
| g__Ureaplasma | 6.432013e-03 | 0.0138012183 |
| g__Lactococcus | 6.366441e-03 | 0.0342174908 |
| g__Lacrimispora | 6.064761e-03 | 0.0098190433 |
| g__Parabacteroides_B | 5.994680e-03 | 0.0100174097 |
| g__CALXRO01 | 5.765728e-03 | 0.0308524398 |
| g__Citrobacter | 5.717035e-03 | 0.0334547978 |
| g__NSJ-61 | 5.560111e-03 | 0.0198839325 |
| g__Breznakia | 5.504528e-03 | 0.0237016548 |
| g__Clostridium_AQ | 5.326190e-03 | 0.0121695008 |
| g__Salmonella | 5.133583e-03 | 0.0184556360 |
| g__Bilophila | 4.983840e-03 | 0.0088456017 |
| g__Hungatella_A | 4.811802e-03 | 0.0095549574 |
| g__MGBC136627 | 4.364796e-03 | 0.0163003355 |
| g__Escherichia | 4.188365e-03 | 0.0266100578 |
| g__UMGS1251 | 4.159842e-03 | 0.0072716746 |
| g__Hungatella | 4.116388e-03 | 0.0190788279 |
| g__Clostridium_Q | 4.011996e-03 | 0.0052127439 |
| g__Brevinema | 4.009680e-03 | 0.0123080282 |
| g__Thomasclavelia | 3.902580e-03 | 0.0109042017 |
| g__Scatousia | 3.649941e-03 | 0.0102727788 |
| g__Enterocloster | 3.620384e-03 | 0.0047339282 |
| g__Mailhella | 3.612079e-03 | 0.0102470769 |
| g__Copromonas | 3.599368e-03 | 0.0050180473 |
| g__Ventrimonas | 3.513355e-03 | 0.0071233140 |
| g__Caccovivens | 3.343316e-03 | 0.0122194941 |
| g__Lawsonia | 3.289015e-03 | 0.0117180962 |
| g__Fournierella | 3.217226e-03 | 0.0062309921 |
| g__Limenecus | 3.163279e-03 | 0.0065731820 |
| g__MGBC133411 | 3.042418e-03 | 0.0074576633 |
| g__Mucinivorans | 2.900095e-03 | 0.0373193950 |
| g__Sarcina | 2.842205e-03 | 0.0171339212 |
| g__Acetatifactor | 2.738138e-03 | 0.0058449483 |
| g__Eisenbergiella | 2.697412e-03 | 0.0068331645 |
| g__Bacteroides_G | 2.682722e-03 | 0.0346150378 |
| g__CAJLXD01 | 2.633818e-03 | 0.0087495783 |
| g__Blautia | 2.558708e-03 | 0.0061407647 |
| g__C-19 | 2.275515e-03 | 0.0048691968 |
| g__Butyricimonas | 2.217402e-03 | 0.0050081070 |
| g__Velocimicrobium | 2.201224e-03 | 0.0066764022 |
| g__CAZU01 | 2.196385e-03 | 0.0065944622 |
| g__MGBC116941 | 2.160169e-03 | 0.0093827415 |
| g__Intestinimonas | 2.081898e-03 | 0.0039382613 |
| g__Negativibacillus | 2.069077e-03 | 0.0055137501 |
| g__Rikenella | 1.985389e-03 | 0.0037067264 |
| g__Phascolarctobacterium | 1.943804e-03 | 0.0050321931 |
| g__RGIG6463 | 1.789843e-03 | 0.0039682804 |
| g__JALFVM01 | 1.742360e-03 | 0.0038623852 |
| g__Oscillibacter | 1.510459e-03 | 0.0024992661 |
| g__Pseudoflavonifractor | 1.502976e-03 | 0.0027706264 |
| g__Acinetobacter | 1.485415e-03 | 0.0097446552 |
| g__Citrobacter_A | 1.392923e-03 | 0.0060348874 |
| g__Staphylococcus | 1.361709e-03 | 0.0050863070 |
| g__RGIG4733 | 1.303790e-03 | 0.0040480788 |
| g__UBA1436 | 1.214079e-03 | 0.0065017516 |
| g__Lachnotalea | 1.208197e-03 | 0.0049169585 |
| g__Ruthenibacterium | 1.202021e-03 | 0.0053632179 |
| g__14-2 | 1.184643e-03 | 0.0096299356 |
| g__Beduini | 1.173950e-03 | 0.0025142208 |
| g__Scatocola | 1.121193e-03 | 0.0044975564 |
| g__Faecisoma | 1.085585e-03 | 0.0055596648 |
| g__Enterococcus_A | 1.083991e-03 | 0.0099150522 |
| g__Faecimonas | 9.880526e-04 | 0.0083079244 |
| g__RGIG9287 | 9.638826e-04 | 0.0093228400 |
| g__CAG-345 | 9.490865e-04 | 0.0060198775 |
| g__Blautia_A | 9.208028e-04 | 0.0029082304 |
| g__RGIG1896 | 8.343610e-04 | 0.0062772562 |
| g__CAG-269 | 7.982448e-04 | 0.0047176841 |
| g__Marseille-P3106 | 7.938268e-04 | 0.0017618833 |
| g__WRHT01 | 6.429563e-04 | 0.0026979819 |
| g__Eggerthella | 6.407882e-04 | 0.0021266467 |
| g__IOR16 | 6.398942e-04 | 0.0020651222 |
| g__Anaerotruncus | 6.265445e-04 | 0.0016948517 |
| g__RUG14156 | 6.151780e-04 | 0.0022147136 |
| g__CHH4-2 | 6.145042e-04 | 0.0019997615 |
| g__Corynebacterium | 5.949530e-04 | 0.0060400054 |
| g__Serratia_A | 5.860616e-04 | 0.0076188009 |
| g__CAG-56 | 4.915432e-04 | 0.0016341973 |
| g__Merdimorpha | 4.736030e-04 | 0.0013057353 |
| g__MGBC140009 | 4.679334e-04 | 0.0024067320 |
| g__CALURL01 | 4.635287e-04 | 0.0016737486 |
| g__Aliarcobacter | 4.604457e-04 | 0.0050324221 |
| g__Scatenecus | 4.520215e-04 | 0.0019760876 |
| g__RGIG8482 | 4.399064e-04 | 0.0029753981 |
| g__Enterobacter | 4.073437e-04 | 0.0041317730 |
| g__Klebsiella | 4.054439e-04 | 0.0048910857 |
| g__Caccenecus | 3.941198e-04 | 0.0017802371 |
| g__Alcaligenes | 3.699430e-04 | 0.0048092594 |
| g__Plesiomonas | 3.633249e-04 | 0.0027105056 |
| g__HGM05232 | 3.617894e-04 | 0.0009776409 |
| g__JAHHSE01 | 3.592120e-04 | 0.0014900895 |
| g__Egerieousia | 3.419549e-04 | 0.0011583178 |
| g__Emergencia | 3.413630e-04 | 0.0017450341 |
| g__Enterococcus_B | 3.352316e-04 | 0.0022266910 |
| g__Stoquefichus | 3.026072e-04 | 0.0020503969 |
| g__Hepatoplasma | 2.989107e-04 | 0.0038858387 |
| g__Paracoccus | 2.959092e-04 | 0.0038468195 |
| g__Moheibacter | 2.771627e-04 | 0.0031476408 |
| g__Scatomorpha | 2.641015e-04 | 0.0010184339 |
| g__UBA7185 | 2.434328e-04 | 0.0014558191 |
| g__Eubacterium | 1.646822e-04 | 0.0006729067 |
| g__Sphingobacterium | 1.505775e-04 | 0.0012460017 |
| g__Devosia | 1.489995e-04 | 0.0015094318 |
| g__Anaerosporobacter | 1.454112e-04 | 0.0012747262 |
| g__Caccomorpha | 1.383122e-04 | 0.0010540604 |
| g__UBA2658 | 1.307451e-04 | 0.0007204965 |
| g__Protoclostridium | 1.277418e-04 | 0.0007646754 |
| g__Angelakisella | 1.268479e-04 | 0.0009221501 |
| g__Cetobacterium_A | 9.310217e-05 | 0.0008718575 |
| g__Rahnella | 6.470705e-05 | 0.0008411917 |
| g__Peptostreptococcus | 2.287338e-05 | 0.0002973539 |
genus_arrange <- genus_summary %>%
group_by(genus) %>%
summarise(mean=sum(relabun)) %>%
filter(genus != "g__")%>%
arrange(-mean) %>%
select(genus) %>%
mutate(genus= sub("^g__", "", genus)) %>%
pull()
genus_summary %>%
left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
mutate(genus= sub("^g__", "", genus)) %>%
filter(genus %in% genus_arrange[1:20]) %>%
mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")
4.3 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))4.3.1 Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric ) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.2 Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.3 Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.4 Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.5 Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.6 Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.7 Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.4 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)4.5 Permanovas
4.5.1 1. Are the wild populations similar?
4.5.1.1 Wild: P.muralis vs P.liolepis
wild <- meta %>%
filter(time_point == "0_Wild")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))
wild_nmds <- sample_metadata %>%
filter(time_point == "0_Wild")4.5.1.3 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.631953 | 0.207198 | 6.795072 | 0.001 |
| Residual | 26 | 6.244347 | 0.792802 | NA | NA |
| Total | 27 | 7.876300 | 1.000000 | NA | NA |
4.5.1.4 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 2.018797 | 0.2566342 | 8.976049 | 0.001 |
| Residual | 26 | 5.847641 | 0.7433658 | NA | NA |
| Total | 27 | 7.866438 | 1.0000000 | NA | NA |
4.5.1.5 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.3786052 | 0.2108638 | 6.947419 | 0.001 |
| Residual | 26 | 1.4168908 | 0.7891362 | NA | NA |
| Total | 27 | 1.7954960 | 1.0000000 | NA | NA |
4.5.1.6 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0787916 | 0.1594096 | 4.930642 | 0.068 |
| Residual | 26 | 0.4154800 | 0.8405904 | NA | NA |
| Total | 27 | 0.4942716 | 1.0000000 | NA | NA |
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
4.5.2 2. Do the antibiotics work?
4.5.2.1 Acclimation vs antibiotics
treat <- meta %>% #antibiotics samples giving error when plotting
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))
treat_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")4.5.2.2 Number of samples used
[1] 53beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)4.5.2.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.033529 | 0.0957379 | 6.818959 | 0.001 |
| species | 1 | 2.327685 | 0.1095866 | 7.805342 | 0.001 |
| individual | 26 | 9.722175 | 0.4577167 | 1.253885 | 0.002 |
| Residual | 24 | 7.157208 | 0.3369589 | NA | NA |
| Total | 52 | 21.240598 | 1.0000000 | NA | NA |
4.5.2.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.168854 | 0.1067496 | 9.686182 | 0.001 |
| species | 1 | 3.090152 | 0.1520953 | 13.800733 | 0.001 |
| individual | 26 | 9.684304 | 0.4766554 | 1.663479 | 0.001 |
| Residual | 24 | 5.373891 | 0.2644996 | NA | NA |
| Total | 52 | 20.317201 | 1.0000000 | NA | NA |
4.5.2.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.0671823 | 0.2208963 | 21.358337 | 0.001 |
| species | 1 | 0.8731487 | 0.0933035 | 9.021460 | 0.001 |
| individual | 26 | 4.0949687 | 0.4375828 | 1.627293 | 0.008 |
| Residual | 24 | 2.3228577 | 0.2482174 | NA | NA |
| Total | 52 | 9.3581574 | 1.0000000 | NA | NA |
4.5.2.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.1253769 | 0.4009905 | 40.2325710 | 0.001 |
| species | 1 | 0.0070451 | 0.0013292 | 0.1333606 | 0.770 |
| individual | 26 | 1.9000410 | 0.3584768 | 1.3833480 | 0.181 |
| Residual | 24 | 1.2678545 | 0.2392035 | NA | NA |
| Total | 52 | 5.3003175 | 1.0000000 | NA | NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
4.5.3 3. Do the FMT work?
4.5.3.1 Comparison between FMT2 vs Post-FMT2
transplant3 <- meta %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")
transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))
transplant3_nmds <- sample_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")4.5.3.2 Number of samples used
[1] 53beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)4.5.3.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9809525 | 0.0778442 | 5.775712 | 0.001 |
| time_point | 1 | 0.7092107 | 0.0562799 | 4.175734 | 0.001 |
| type | 1 | 1.4479860 | 0.1149060 | 8.525540 | 0.001 |
| individual | 25 | 6.9157236 | 0.5488022 | 1.628753 | 0.001 |
| Residual | 15 | 2.5476146 | 0.2021678 | NA | NA |
| Total | 43 | 12.6014875 | 1.0000000 | NA | NA |
4.5.3.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.1101711 | 0.0895457 | 7.879642 | 0.001 |
| time_point | 1 | 0.7363935 | 0.0593971 | 5.226687 | 0.001 |
| type | 1 | 1.9027421 | 0.1534740 | 13.505059 | 0.001 |
| individual | 25 | 6.5351386 | 0.5271203 | 1.855374 | 0.001 |
| Residual | 15 | 2.1133659 | 0.1704628 | NA | NA |
| Total | 43 | 12.3978112 | 1.0000000 | NA | NA |
4.5.3.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1341431 | 0.1005265 | 6.921553 | 0.001 |
| time_point | 1 | 0.1159136 | 0.0868654 | 5.980943 | 0.001 |
| type | 1 | 0.1423573 | 0.1066822 | 7.345390 | 0.001 |
| individual | 25 | 0.6512838 | 0.4880704 | 1.344205 | 0.079 |
| Residual | 15 | 0.2907075 | 0.2178554 | NA | NA |
| Total | 43 | 1.3344053 | 1.0000000 | NA | NA |
4.5.3.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0031096 | -0.0030953 | -0.1177247 | 0.818 |
| time_point | 1 | -0.0045864 | -0.0045653 | -0.1736359 | 0.827 |
| type | 1 | 0.0775554 | 0.0771987 | 2.9361453 | 0.143 |
| individual | 25 | 0.5385511 | 0.5360739 | 0.8155530 | 0.631 |
| Residual | 15 | 0.3962105 | 0.3943880 | NA | NA |
| Total | 43 | 1.0046211 | 1.0000000 | NA | NA |
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylo_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_func_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")
4.5.4 4. Are there differences between the control and the treatment group?
4.5.4.2 Number of samples used
[1] 27beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)4.5.4.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.6488906 | 0.0757032 | 2.115638 | 0.002 |
| type | 1 | 0.5615418 | 0.0655126 | 1.830847 | 0.006 |
| Residual | 24 | 7.3610760 | 0.8587842 | NA | NA |
| Total | 26 | 8.5715084 | 1.0000000 | NA | NA |
4.5.4.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.7984743 | 0.104299 | 3.065171 | 0.001 |
| type | 1 | 0.6051778 | 0.079050 | 2.323148 | 0.004 |
| Residual | 24 | 6.2519772 | 0.816651 | NA | NA |
| Total | 26 | 7.6556293 | 1.000000 | NA | NA |
4.5.4.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0700374 | 0.0635900 | 1.6594454 | 0.140 |
| type | 1 | 0.0184254 | 0.0167292 | 0.4365646 | 0.781 |
| Residual | 24 | 1.0129278 | 0.9196808 | NA | NA |
| Total | 26 | 1.1013906 | 1.0000000 | NA | NA |
4.5.4.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0048931 | -0.0058792 | -0.1628739 | 0.837 |
| type | 1 | 0.1161466 | 0.1395538 | 3.8660898 | 0.112 |
| Residual | 24 | 0.7210172 | 0.8663254 | NA | NA |
| Total | 26 | 0.8322707 | 1.0000000 | NA | NA |
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")
4.5.4.4 Number of samples used
[1] 28beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)4.5.4.4.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.8966107 | 0.1132233 | 3.515588 | 0.001 |
| type | 1 | 0.6463814 | 0.0816245 | 2.534445 | 0.002 |
| Residual | 25 | 6.3759661 | 0.8051521 | NA | NA |
| Total | 27 | 7.9189582 | 1.0000000 | NA | NA |
4.5.4.4.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9398968 | 0.1224370 | 4.112304 | 0.001 |
| type | 1 | 1.0227481 | 0.1332297 | 4.474800 | 0.001 |
| Residual | 25 | 5.7139316 | 0.7443333 | NA | NA |
| Total | 27 | 7.6765766 | 1.0000000 | NA | NA |
4.5.4.4.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1200820 | 0.1454923 | 4.647187 | 0.001 |
| type | 1 | 0.0592745 | 0.0718175 | 2.293931 | 0.027 |
| Residual | 25 | 0.6459931 | 0.7826902 | NA | NA |
| Total | 27 | 0.8253496 | 1.0000000 | NA | NA |
4.5.4.4.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0116738 | 0.0168395 | 0.4229082 | 0.551 |
| type | 1 | -0.0085273 | -0.0123007 | -0.3089208 | 0.909 |
| Residual | 25 | 0.6900904 | 0.9954612 | NA | NA |
| Total | 27 | 0.6932368 | 1.0000000 | NA | NA |
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")